/****************************************************************************
*                                                                           *
* scatter.cpp - scattering simulator			                            *
*                                                                           *
* ALaRI 2008												                *
*                                                                           *
****************************************************************************/

#include "scatter.h"	// - me
#include "matrepository.h"
#include "sensorrepository.h"
#include "sensor.h"

#include <memory>

#include <QPoint> 
#include <QRect> 
#include <QColor>
#include <QFuture>
#include <QtConcurrentRun>

/****************************************************************************
* CScatSimulator - scattering parflow simulator				                *
****************************************************************************/

/****************************************************************************
* Construction/Destruction													*
****************************************************************************/
CScatSimulator::CScatSimulator(const int& nWidth, const int& nHeight, const QImage& qImageFloor) :
	m_pSS(NULL),
	m_pB1(NULL),
	m_pB2(NULL),
	m_pScatMatrix(NULL),
	m_nWidth(nWidth),
	m_nHeight(nHeight),
	m_qImageFloor(qImageFloor)
{}

CScatSimulator::~CScatSimulator()
{
	// destroy floorplan
	freeSimulationData();
}

/****************************************************************************
* Public																	*
****************************************************************************/
float CScatSimulator::getFieldAt(const QPoint& point, CMaterial *pMaterial) const
{
	// simulation should have finished
	assert(m_pSS != NULL);
	assert(NULL != pMaterial);
	if (pMaterial == NULL)
	{
		/****************************************************************************
		* TODO: throw an exception here !!!!! 										*
		****************************************************************************/
		return -1.0;
	}

	// valid point
	assert(point.x() >= 0 && point.y() >= 0);
	assert(point.x() <= m_nWidth && point.y() <= m_nHeight);

	// get point in simulation space
	SSimPoint& simPoint = m_pSS[point.y() * m_nWidth + point.x()];

	return computeField(simPoint, pMaterial);

}

float CScatSimulator::computeField(SSimPoint& simPoint, CMaterial *pMaterial) const
{
	// calculate auxiliary var
	std::complex<float> cfltKN2 = pMaterial->calculateKN2();
//	return abs(cfltKN2 * (simPoint.m_cfltIN + simPoint.m_cfltIE +
//		simPoint.m_cfltIS + simPoint.m_cfltIW));
	return pow(abs(cfltKN2 * (simPoint.m_cfltIN + simPoint.m_cfltIE +
		simPoint.m_cfltIS + simPoint.m_cfltIW)), 2.f);
}

void CScatSimulator::elaborateFloorPlan()
{
	if (m_pScatMatrix == NULL)
	{
		// create rectangle for entire area
		QRect rect(QPoint(0,0), QPoint(m_nWidth, m_nHeight));
		// get scattering matrix blocks - will cause elaboration of entire space
		m_pScatMatrix = getScatMatrix(rect);
		m_pScatMatrix->getBlock(CScatMatrix::eBlock_ILT, true);
		if (true == m_bCancelSimulation)
			return;
		m_pScatMatrix->getBlock(CScatMatrix::eBlock_IRB, true);
	}
}

void CScatSimulator::putSensor(CSensor *pSensor) 
{
	int nPosX = pSensor->x();
	int nPosY = pSensor->y();

	assert(nPosX >= 0 && nPosX < m_nWidth);
	assert(nPosY >= 0 && nPosY < m_nHeight);

	freeSimSpace();

	// calculate initial value of the field
	if (m_pSS == NULL)
		allocateSimSpace();
	// clean it
	for (int nIndex = m_nWidth * m_nHeight; --nIndex >= 0; )
		m_pSS[nIndex].m_cfltIN = m_pSS[nIndex].m_cfltIE =
		m_pSS[nIndex].m_cfltIS = m_pSS[nIndex].m_cfltIW = 0.f;

	assert(nPosX >= 0 && nPosY >= 0);
	assert(nPosX < m_nWidth && nPosY < m_nHeight);

	m_pSS[(nPosY - 1) * m_nWidth + nPosX].m_cfltIS = 
	m_pSS[nPosY * m_nWidth + (nPosX - 1)].m_cfltIE =
	m_pSS[nPosY * m_nWidth + (nPosX + 1)].m_cfltIW =
	m_pSS[(nPosY + 1) * m_nWidth + nPosX].m_cfltIN = 100.;
	
	QRgb rgbColor = m_qImageFloor.pixel(pSensor->x(), pSensor->y());
	CMaterial *pMaterial = CMaterialRepository::inst()->getMaterial(rgbColor);
	assert(NULL != pMaterial);
	SSimPoint simPoint = {100., 100., 100., 100.};
	pSensor->setRefField(computeField(simPoint, pMaterial)); 
}

void CScatSimulator::runSimulation(bool bElaborate)
{
	m_bCancelSimulation = false;
	m_simWatcher.disconnect();
	connect(&m_simWatcher, SIGNAL(finished()), this,  SIGNAL(simulationDone()));
	// create concurrent thread which will run the simulation
	QFuture<void> future = QtConcurrent::run(this, &CScatSimulator::simulate, bElaborate);
	m_simWatcher.setFuture(future);
}

void CScatSimulator::simulate(bool bElaborate)
{
	if (bElaborate)
	{
		emit doing(tr("Elaborating Floor Plan..."));  
		freeElaborationData();
		elaborateFloorPlan();
	}
	QList<CSensor*> qListSensors = CSensorRepository::inst()->toQList();
	QListIterator<CSensor*> it(qListSensors);

	while (it.hasNext())
	{
		CSensor *pSensor = it.next();
		assert(NULL != pSensor);
		if (NULL == pSensor)
			// TODO: throw an exception here
			return; 
		if (true == pSensor->isFieldUpdated() && false == bElaborate) 
			continue;
		putSensor(pSensor); 

		if (true == m_bCancelSimulation)
			return;
		
		QRect rect(QPoint(0, 0), QPoint(m_nWidth, m_nHeight));
		emit doing(tr("Simulation Phase 1..."));  
		phase1UpdateInnerFlows(rect, m_pScatMatrix);
		if (true == m_bCancelSimulation)
			return;

		emit doing(tr("Simulation Phase 2..."));  
		phase2UpdateInnerFlows(rect, m_pScatMatrix);
		if (true == m_bCancelSimulation)
			return;
		
		QPoint qPt;
		for (qPt.ry() = 0; qPt.ry() < m_nHeight; qPt.ry()++)
		{
			for (qPt.rx() = 0; qPt.rx() < m_nWidth; qPt.rx()++)
			{
				QRgb rgbColor = m_qImageFloor.pixel(qPt.x(), qPt.y());
				float fltValue = abs(10.*log10(getFieldAt(qPt, 
					CMaterialRepository::inst()->getMaterial(rgbColor))/pSensor->getRefField())); 
				int byValue = (int)(fltValue);
				if (byValue > 255)
					byValue = 255;
				pSensor->setFieldAt(qPt, byValue);
			}
		}
		pSensor->setFieldUpdated(true);
	}
}

const CScatMatrix *CScatSimulator::getScatMatrix(const QRect& rect)
{
	assert(rect.left() >= 0 && rect.top() >= 0);
	assert(rect.right() > rect.left() && rect.bottom() > rect.top());
	assert(rect.right() <= m_nWidth && rect.bottom() <= m_nHeight);

	// level 0 matrix?
	if (rect.width() - 1 == 1 && rect.height() - 1 == 1)
		return getScatMatrixL0(rect.topLeft());

	// get scat matrices from children
	QRect rectChild; const CScatMatrix *pLOrT, *pROrB;
	bool bHorizontal = rect.width() >= rect.height();

	// see if we need to devide vertically or horizontally
	if (bHorizontal)
	{
		// left child
		rectChild.setTopLeft(rect.topLeft());
		rectChild.setRight((rect.left() + rect.right()) / 2);
		rectChild.setBottom(rect.bottom());
		pLOrT = getScatMatrix(rectChild);

		// right child
		rectChild.setLeft(rectChild.right());
		rectChild.setRight(rect.right());
		pROrB = getScatMatrix(rectChild);
	}
	else
	{
		// top child
		rectChild.setTopLeft(rect.topLeft());
		rectChild.setBottom((rect.top() + rect.bottom()) / 2);
		rectChild.setRight(rect.right());
		pLOrT = getScatMatrix(rectChild);

		// bottom child
		rectChild.setTop(rectChild.bottom());
		rectChild.setBottom(rect.bottom());
		pROrB = getScatMatrix(rectChild);
	}

	// now lookup the entry in cache to see whether we already have this matrix
	CScatMatrixLP *pMatrix; 
	__SScatKey key = { pLOrT, pROrB, bHorizontal };
	pMatrix = m_mapScatMatrices[key];
	if (pMatrix == NULL)
	{
		// create the new one
		std::auto_ptr<CScatMatrixLP> spMatrix(new CScatMatrixLP(pLOrT, pROrB, bHorizontal));
		// add it to the cache
		m_mapScatMatrices[key] = spMatrix.get();
		// and detach
		pMatrix = spMatrix.release();
	}

	return pMatrix;
}


const CScatMatrix* CScatSimulator::getScatMatrixL0(const QPoint& point)
{
	assert(point.x() >= 0 && point.y() >= 0);
	assert(point.x() <= m_nWidth && point.y() <= m_nHeight);

	QRgb rgbColor = m_qImageFloor.pixel(point.x(), point.y());
	CMaterial *pMaterial = CMaterialRepository::inst()->getMaterial(rgbColor);
	
	if (pMaterial == NULL)
	{
		/****************************************************************************
		* throw an exception here !!!!! 											*
		****************************************************************************/
		assert(NULL);
		return NULL;
	}

	CScatMatrix *pScatMatrix = pMaterial->getScatMatrix();
	if (pScatMatrix == NULL)
	{
		/****************************************************************************
		* throw an exception here !!!!! 											*
		****************************************************************************/
		return NULL;
	}
	
	return pScatMatrix;
}

void CScatSimulator::cancelSimulation()
{
	m_bCancelSimulation = true;
	m_simWatcher.waitForFinished();
	m_simWatcher.disconnect();
	freeSimulationData();
}

void CScatSimulator::freeSimulationData()
{
	// destroy simulation space
	freeSimSpace();
	// free elaboration data
	freeElaborationData();
}

void CScatSimulator::freeElaborationData()
{
	if (m_pScatMatrix != NULL)
	{
		// delete scattering matrices for the floorplan
		std::hash_map<__SScatKey, CScatMatrixLP*, 
			__SKeyHasher, __SKeysComparator>::iterator itMatrices = m_mapScatMatrices.begin();

		while(itMatrices != m_mapScatMatrices.end())
		{
			SafeFree(itMatrices->second);
			itMatrices++;
		}
		
		m_mapScatMatrices.clear();
		m_pScatMatrix = NULL;
	}
}

void CScatSimulator::allocateSimSpace()
{
	// should not be allocated
	if (m_pSS == NULL)
	{
		assert(m_nWidth > 0 && m_nHeight > 0);

		// allocate space for simulation
		std::auto_ptr<SSimPoint> spSS(new SSimPoint[m_nWidth * m_nHeight]);

		int nSize; // allocate space for largest update array
		if (m_nWidth > m_nHeight)
			nSize = (m_nWidth / 2 + 1 > m_nHeight) ? (m_nWidth / 2 + 1) : m_nHeight;
		else
			nSize = (m_nHeight / 2 + 1 > m_nWidth) ? (m_nHeight / 2 + 1) : m_nWidth;

		std::auto_ptr<std::complex<float> > spB1(new std::complex<float>[nSize]);
		std::auto_ptr<std::complex<float> > spB2(new std::complex<float>[nSize]);

		// save them
		m_pSS = spSS.release();
		m_pB1 = spB1.release();
		m_pB2 = spB2.release();
	}
}

void CScatSimulator::freeSimSpace()
{
	// delete and reset if there
	if (m_pSS != NULL)
	{
		SafeFree(m_pSS);
		SafeFree(m_pB1);
		SafeFree(m_pB2);
	}
}

void CScatSimulator::phase1UpdateInnerFlows(const QRect& rect, const CScatMatrix *pMatrix)
{
	assert(rect.left() >= 0 && rect.top() >= 0);
	assert(rect.right() > rect.left() && rect.bottom() > rect.top());
	assert(rect.width() - 1 + rect.height() - 1 > 2);
	assert(rect.right() <= m_nWidth && rect.bottom() <= m_nHeight);
	assert(pMatrix != NULL);

	// get all 4 blocks of the I Matrix
	const CScatBlock *pILT = pMatrix->getBlock(CScatMatrix::eBlock_ILT, true);
	const CScatBlock *pIRT = pMatrix->getBlock(CScatMatrix::eBlock_IRT, true);
	const CScatBlock *pILB = pMatrix->getBlock(CScatMatrix::eBlock_ILB, true);
	const CScatBlock *pIRB = pMatrix->getBlock(CScatMatrix::eBlock_IRB, true);

	int nRow, nColumn;
	QRect rectChild; bool bHorizontal = rect.width() >= rect.height();

	// see if we need to devide vertically or horizontally
	if (bHorizontal)
	{
		assert(pILT->getHeight() == rect.height() - 1);

		// left child
		rectChild.setTopLeft(rect.topLeft());
		rectChild.setRight((rect.left() + rect.right()) / 2);
		rectChild.setBottom(rect.bottom());
		if (rectChild.width() - 1 != 1 || rectChild.height() - 1 != 1)
			phase1UpdateEquivalentSource(rectChild, pMatrix->getLOrT());

		// right child
		rectChild.setLeft(rectChild.right());
		rectChild.setRight(rect.right());
		if (rectChild.width() - 1 != 1 || rectChild.height() - 1 != 1)
			phase1UpdateEquivalentSource(rectChild, pMatrix->getROrB());

		// update accordingly
		for (nRow = pILT->getHeight(); --nRow >= 0; )
			for (m_pB1[nRow] = 0.f, m_pB2[nRow] = 0.f, nColumn = pILT->getWidth(); --nColumn >= 0;)
			{
				// I
				m_pB1[nRow] += pILT->valueAt(nRow, nColumn) *
						m_pSS[(rectChild.top() + nColumn) * m_nWidth + rectChild.left()].m_cfltIW +
					pIRT->valueAt(nRow, nColumn) *
						m_pSS[(rectChild.top() + nColumn) * m_nWidth + rectChild.left() - 1].m_cfltIE;
				// J
				m_pB2[nRow] += pILB->valueAt(nRow, nColumn) *
						m_pSS[(rectChild.top() + nColumn) * m_nWidth + rectChild.left()].m_cfltIW +
					pIRB->valueAt(nRow, nColumn) *
						m_pSS[(rectChild.top() + nColumn) * m_nWidth + rectChild.left() - 1].m_cfltIE;
			}

		// write the results back
		for (nRow = pILT->getHeight(); --nRow >= 0; )
			m_pSS[(rectChild.top() + nRow) * m_nWidth + rectChild.left() - 1].m_cfltIE = m_pB1[nRow],
			m_pSS[(rectChild.top() + nRow) * m_nWidth + rectChild.left()].m_cfltIW = m_pB2[nRow];
	}
	else
	{
		assert(pILT->getHeight() == rect.width() - 1);

		// top child
		rectChild.setTopLeft(rect.topLeft());
		rectChild.setBottom((rect.top() + rect.bottom()) / 2);
		rectChild.setRight(rect.right());
		if (rectChild.width() - 1 != 1 || rectChild.height() - 1 != 1)
			phase1UpdateEquivalentSource(rectChild, pMatrix->getLOrT());

		// bottom child
		rectChild.setTop(rectChild.bottom());
		rectChild.setBottom(rect.bottom());
		if (rectChild.width() - 1 != 1 || rectChild.height() - 1 != 1)
			phase1UpdateEquivalentSource(rectChild, pMatrix->getROrB());

		// update accordingly
		for (nRow = pILT->getHeight(); --nRow >= 0; )
			for (m_pB1[nRow] = 0.f, m_pB2[nRow] = 0.f, nColumn = pILT->getWidth(); --nColumn >= 0;)
			{
				// I
				m_pB1[nRow] += pILT->valueAt(nRow, nColumn) *
						m_pSS[rectChild.top() * m_nWidth + rectChild.left() + nColumn].m_cfltIN +
					pIRT->valueAt(nRow, nColumn) *
						m_pSS[(rectChild.top() - 1) * m_nWidth + rectChild.left() + nColumn].m_cfltIS; 
				// J
				m_pB2[nRow] += pILB->valueAt(nRow, nColumn) *
						m_pSS[rectChild.top() * m_nWidth + rectChild.left() + nColumn].m_cfltIN +
					pIRB->valueAt(nRow, nColumn) *
						m_pSS[(rectChild.top() - 1) * m_nWidth + rectChild.left() + nColumn].m_cfltIS; 
			}

		// write the results back
		for (nRow = pILT->getHeight(); --nRow >= 0; )
			m_pSS[(rectChild.top() - 1) * m_nWidth + rectChild.left() + nRow].m_cfltIS = m_pB1[nRow],
			m_pSS[rectChild.top() * m_nWidth + rectChild.left() + nRow].m_cfltIN = m_pB2[nRow];
	}
}

void CScatSimulator::phase1UpdateEquivalentSource(const QRect& rect, const CScatMatrix *pMatrix)
{
	assert(rect.left() >= 0 && rect.top() >= 0);
	assert(rect.right() > rect.left() && rect.bottom() > rect.top());
	assert(rect.width() - 1 + rect.height() - 1 > 2);
	assert(rect.right() <= m_nWidth && rect.bottom() <= m_nHeight);
	assert(pMatrix != NULL);

	// first we need to update internal flows
	phase1UpdateInnerFlows(rect, pMatrix);

	// get child matrixes
	const CScatMatrix *pLOrT = pMatrix->getLOrT();
	const CScatMatrix *pROrB = pMatrix->getROrB();

	// other vars
	int nRow, nColumn; const CScatBlock* pB; 	
	QPoint centerPt(rect.center());

	// see if we need to devide vertically or horizontally
	bool bHorizontal = rect.width() >= rect.height();
	if (bHorizontal)
	{
		// IEN
		if (rect.top() > 0)
		{
			pB = pLOrT->getBlock(CScatMatrix::eBlock_E2N, true);
			for (nRow = pB->getHeight(); --nRow >= 0;)
				for (nColumn = pB->getWidth(); --nColumn >= 0; )
					m_pSS[(rect.top() - 1) * m_nWidth + rect.left() + nRow].m_cfltIS +=
						pB->valueAt(nRow, nColumn) *
							m_pSS[(rect.top() + nColumn) * m_nWidth + centerPt.x() - 1].m_cfltIE;
		}
		// IES
		if (rect.bottom() < m_nHeight)
		{
			pB = pLOrT->getBlock(CScatMatrix::eBlock_E2S, true);
			for (nRow = pB->getHeight(); --nRow >= 0;)
				for (nColumn = pB->getWidth(); --nColumn >= 0; )
					m_pSS[rect.bottom() * m_nWidth + rect.left() + nRow].m_cfltIN +=
						pB->valueAt(nRow, nColumn) *
							m_pSS[(rect.top() + nColumn) * m_nWidth + centerPt.x() - 1].m_cfltIE;
		}
		// IEW
		if (rect.left() > 0)
		{
			pB = pLOrT->getBlock(CScatMatrix::eBlock_E2W, true);
			for (nRow = pB->getHeight(); --nRow >= 0;)
				for (nColumn = pB->getWidth(); --nColumn >= 0; )
					m_pSS[(rect.top() + nRow) * m_nWidth + rect.left() - 1].m_cfltIE +=
						pB->valueAt(nRow, nColumn) *
							m_pSS[(rect.top() + nColumn) * m_nWidth + centerPt.x() - 1].m_cfltIE;
		}
		// JWN
		if (rect.top() > 0)
		{
			pB = pROrB->getBlock(CScatMatrix::eBlock_W2N, true);
			for (nRow = pB->getHeight(); --nRow >= 0;)
				for (nColumn = pB->getWidth(); --nColumn >= 0; )
					m_pSS[(rect.top() - 1) * m_nWidth + centerPt.x() + nRow].m_cfltIS +=
						pB->valueAt(nRow, nColumn) *
							m_pSS[(rect.top() + nColumn) * m_nWidth + centerPt.x()].m_cfltIW;
		}
		// JWE
		if (rect.right() < m_nWidth)
		{
			pB = pROrB->getBlock(CScatMatrix::eBlock_W2E, true);
			for (nRow = pB->getHeight(); --nRow >= 0;)
				for (nColumn = pB->getWidth(); --nColumn >= 0; )
					m_pSS[(rect.top() + nRow) * m_nWidth + rect.right()].m_cfltIW +=
						pB->valueAt(nRow, nColumn) *
							m_pSS[(rect.top() + nColumn) * m_nWidth + centerPt.x()].m_cfltIW;
		}
		// JWS
		if (rect.bottom() < m_nHeight)
		{
			pB = pROrB->getBlock(CScatMatrix::eBlock_W2S, true);
			for (nRow = pB->getHeight(); --nRow >= 0;)
				for (nColumn = pB->getWidth(); --nColumn >= 0; )
					m_pSS[rect.bottom() * m_nWidth + centerPt.x() + nRow].m_cfltIN +=
						pB->valueAt(nRow, nColumn) *
							m_pSS[(rect.top() + nColumn) * m_nWidth + centerPt.x()].m_cfltIW;
		}
	}
	else // vertical
	{
		// ISN
		if (rect.top() > 0)
		{
			pB = pLOrT->getBlock(CScatMatrix::eBlock_S2N, true);
			for (nRow = pB->getHeight(); --nRow >= 0;)
				for (nColumn = pB->getWidth(); --nColumn >= 0; )
					m_pSS[(rect.top() - 1) * m_nWidth + rect.left() + nRow].m_cfltIS +=
						pB->valueAt(nRow, nColumn) *
							m_pSS[(centerPt.y() - 1) * m_nWidth + rect.left() + nColumn].m_cfltIS;
		}
		// ISE
		if (rect.right() < m_nWidth)
		{
			pB = pLOrT->getBlock(CScatMatrix::eBlock_S2E, true);
			for (nRow = pB->getHeight(); --nRow >= 0;)
				for (nColumn = pB->getWidth(); --nColumn >= 0; )
					m_pSS[(rect.top() + nRow) * m_nWidth + rect.right()].m_cfltIW +=
						pB->valueAt(nRow, nColumn) *
							m_pSS[(centerPt.y() - 1) * m_nWidth + rect.left() + nColumn].m_cfltIS;
		}
		// ISW
		if (rect.left() > 0)
		{
			pB = pLOrT->getBlock(CScatMatrix::eBlock_S2W, true);
			for (nRow = pB->getHeight(); --nRow >= 0;)
				for (nColumn = pB->getWidth(); --nColumn >= 0; )
					m_pSS[(rect.top() + nRow) * m_nWidth + rect.left() - 1].m_cfltIE +=
						pB->valueAt(nRow, nColumn) *
							m_pSS[(centerPt.y() - 1) * m_nWidth + rect.left() + nColumn].m_cfltIS;
		}
		// JNE
		if (rect.right() < m_nWidth)
		{
			pB = pROrB->getBlock(CScatMatrix::eBlock_N2E, true);
			for (nRow = pB->getHeight(); --nRow >= 0;)
				for (nColumn = pB->getWidth(); --nColumn >= 0; )
					m_pSS[(centerPt.y() + nRow) * m_nWidth + rect.right()].m_cfltIW +=
						pB->valueAt(nRow, nColumn) *
							m_pSS[centerPt.y() * m_nWidth + rect.left() + nColumn].m_cfltIN;
		}
		// JNS
		if (rect.bottom() < m_nHeight)
		{
			pB = pROrB->getBlock(CScatMatrix::eBlock_N2S, true);
			for (nRow = pB->getHeight(); --nRow >= 0;)
				for (nColumn = pB->getWidth(); --nColumn >= 0; )
					m_pSS[rect.bottom() * m_nWidth + rect.left() + nRow].m_cfltIN +=
						pB->valueAt(nRow, nColumn) *
							m_pSS[centerPt.y() * m_nWidth + rect.left() + nColumn].m_cfltIN;
		}
		// JNW
		if (rect.left() > 0)
		{
			pB = pROrB->getBlock(CScatMatrix::eBlock_N2W, true);
			for (nRow = pB->getHeight(); --nRow >= 0;)
				for (nColumn = pB->getWidth(); --nColumn >= 0; )
					m_pSS[(centerPt.y() + nRow) * m_nWidth + rect.left() - 1].m_cfltIE +=
						pB->valueAt(nRow, nColumn) *
							m_pSS[centerPt.y() * m_nWidth + rect.left() + nColumn].m_cfltIN;
		}
	}
}

void CScatSimulator::phase2UpdateInnerFlows(const QRect& rect, const CScatMatrix *pMatrix)
{
	assert(rect.left() >= 0 && rect.top() >= 0);
	assert(rect.right() > rect.left() && rect.bottom() > rect.top());
	assert(rect.width() - 1 + rect.height() - 1 > 2);
	assert(rect.right() <= m_nWidth && rect.bottom() <= m_nHeight);
	assert(pMatrix != NULL);

	// get child matrixes
	const CScatMatrix *pLOrT = pMatrix->getLOrT();
	const CScatMatrix *pROrB = pMatrix->getROrB();

	// get all 4 blocks of the I Matrix
	const CScatBlock *pILT = pMatrix->getBlock(CScatMatrix::eBlock_ILT, true);
	const CScatBlock *pIRT = pMatrix->getBlock(CScatMatrix::eBlock_IRT, true);
	const CScatBlock *pILB = pMatrix->getBlock(CScatMatrix::eBlock_ILB, true);
	const CScatBlock *pIRB = pMatrix->getBlock(CScatMatrix::eBlock_IRB, true);

	// other vars
	int nRow, nColumn; const CScatBlock* pB; 	
	QRect rectChild; QPoint centerPt(rect.center());

	// see if we need to devide vertically or horizontally
	bool bHorizontal = rect.width() >= rect.height();
	if (bHorizontal)
	{
		// cleanup buffers first
		for (nRow = rect.height() - 1; --nRow >= 0; )
			m_pB1[nRow] = m_pB2[nRow] = 0.f;

		// INE
		pB = pLOrT->getBlock(CScatMatrix::eBlock_N2E, true);
		for (nRow = pB->getHeight(); --nRow >= 0;)
			for (nColumn = pB->getWidth(); --nColumn >= 0; )
				m_pB1[nRow] += pB->valueAt(nRow, nColumn) *
					m_pSS[rect.top() * m_nWidth + rect.left() + nColumn].m_cfltIN;
		// IWE
		pB = pLOrT->getBlock(CScatMatrix::eBlock_W2E, true);
		for (nRow = pB->getHeight(); --nRow >= 0;)
			for (nColumn = pB->getWidth(); --nColumn >= 0; )
				m_pB1[nRow] += pB->valueAt(nRow, nColumn) *
					m_pSS[(rect.top() + nColumn) * m_nWidth + rect.left()].m_cfltIW;
		// ISE
		pB = pLOrT->getBlock(CScatMatrix::eBlock_S2E, true);
		for (nRow = pB->getHeight(); --nRow >= 0;)
			for (nColumn = pB->getWidth(); --nColumn >= 0; )
				m_pB1[nRow] += pB->valueAt(nRow, nColumn) *
					m_pSS[(rect.bottom() - 1) * m_nWidth + rect.left() + nColumn].m_cfltIS;
		// JNW
		pB = pROrB->getBlock(CScatMatrix::eBlock_N2W, true);
		for (nRow = pB->getHeight(); --nRow >= 0;)
			for (nColumn = pB->getWidth(); --nColumn >= 0; )
				m_pB2[nRow] += pB->valueAt(nRow, nColumn) *
					m_pSS[rect.top() * m_nWidth + centerPt.x() + nColumn].m_cfltIN;
		// JEW
		pB = pROrB->getBlock(CScatMatrix::eBlock_E2W, true);
		for (nRow = pB->getHeight(); --nRow >= 0;)
			for (nColumn = pB->getWidth(); --nColumn >= 0; )
				m_pB2[nRow] += pB->valueAt(nRow, nColumn) *
					m_pSS[(rect.top() + nColumn) * m_nWidth + rect.right() - 1].m_cfltIE;
		// JSW
		pB = pROrB->getBlock(CScatMatrix::eBlock_S2W, true);
		for (nRow = pB->getHeight(); --nRow >= 0;)
			for (nColumn = pB->getWidth(); --nColumn >= 0; )
				m_pB2[nRow] += pB->valueAt(nRow, nColumn) *
					m_pSS[(rect.bottom() - 1) * m_nWidth + centerPt.x() + nColumn].m_cfltIS;

		// now new internal source updates are in buffers
		// we have to pass them through I matrix
		// ILT, ILB
		for (nRow = pILT->getHeight(); --nRow >= 0; )
			for (nColumn = pILT->getWidth(); --nColumn >= 0;)
			{
				// I
				m_pSS[(rect.top() + nRow) * m_nWidth + centerPt.x() - 1].m_cfltIE +=
					pILT->valueAt(nRow, nColumn) * m_pB1[nColumn] +
						pIRT->valueAt(nRow, nColumn) * m_pB2[nColumn];
				// J
				m_pSS[(rect.top() + nRow) * m_nWidth + centerPt.x()].m_cfltIW +=
					pILB->valueAt(nRow, nColumn) * m_pB1[nColumn] +
						pIRB->valueAt(nRow, nColumn) * m_pB2[nColumn];
			}

		// no update children - left one
		rectChild.setTopLeft(rect.topLeft());
		rectChild.setRight(centerPt.x());
		rectChild.setBottom(rect.bottom());
		if (rectChild.width() - 1 != 1 || rectChild.height() - 1 != 1)
			phase2UpdateInnerFlows(rectChild, pLOrT);

		// right child
		rectChild.setLeft(rectChild.right());
		rectChild.setRight(rect.right());
		if (rectChild.width() - 1 != 1 || rectChild.height() - 1 != 1)
			phase2UpdateInnerFlows(rectChild, pROrB);
	}
	else // vertical
	{
		// cleanup buffers first
		for (nRow = rect.width() - 1; --nRow >= 0; )
			m_pB1[nRow] = m_pB2[nRow] = 0.f;

		// INS
		pB = pLOrT->getBlock(CScatMatrix::eBlock_N2S, true);
		for (nRow = pB->getHeight(); --nRow >= 0;)
			for (nColumn = pB->getWidth(); --nColumn >= 0; )
				m_pB1[nRow] += pB->valueAt(nRow, nColumn) *
					m_pSS[rect.top() * m_nWidth + rect.left() + nColumn].m_cfltIN;
		// IES
		pB = pLOrT->getBlock(CScatMatrix::eBlock_E2S, true);
		for (nRow = pB->getHeight(); --nRow >= 0;)
			for (nColumn = pB->getWidth(); --nColumn >= 0; )
				m_pB1[nRow] += pB->valueAt(nRow, nColumn) *
					m_pSS[(rect.top() + nColumn) * m_nWidth + rect.right() - 1].m_cfltIE;
		// IWS
		pB = pLOrT->getBlock(CScatMatrix::eBlock_W2S, true);
		for (nRow = pB->getHeight(); --nRow >= 0;)
			for (nColumn = pB->getWidth(); --nColumn >= 0; )
				m_pB1[nRow] += pB->valueAt(nRow, nColumn) *
					m_pSS[(rect.top() + nColumn) * m_nWidth + rect.left()].m_cfltIW;
		// JEN
		pB = pROrB->getBlock(CScatMatrix::eBlock_E2N, true);
		for (nRow = pB->getHeight(); --nRow >= 0;)
			for (nColumn = pB->getWidth(); --nColumn >= 0; )
				m_pB2[nRow] += pB->valueAt(nRow, nColumn) *
					m_pSS[(centerPt.y() + nColumn) * m_nWidth + rect.right() - 1].m_cfltIE;
		// JSN
		pB = pROrB->getBlock(CScatMatrix::eBlock_S2N, true);
		for (nRow = pB->getHeight(); --nRow >= 0;)
			for (nColumn = pB->getWidth(); --nColumn >= 0; )
				m_pB2[nRow] += pB->valueAt(nRow, nColumn) *
					m_pSS[(rect.bottom() - 1) * m_nWidth + rect.left() + nColumn].m_cfltIS;
		// JWN
		pB = pROrB->getBlock(CScatMatrix::eBlock_W2N, true);
		for (nRow = pB->getHeight(); --nRow >= 0;)
			for (nColumn = pB->getWidth(); --nColumn >= 0; )
				m_pB2[nRow] += pB->valueAt(nRow, nColumn) *
					m_pSS[(centerPt.y() + nColumn) * m_nWidth + rect.left()].m_cfltIW;

		// ILT, ILB
		for (nRow = pILT->getHeight(); --nRow >= 0; )
			for (nColumn = pILT->getWidth(); --nColumn >= 0;)
			{
				// I - i changed nRow to nColumn
				m_pSS[(centerPt.y() - 1) * m_nWidth + rect.left() + nRow].m_cfltIS +=
					pILT->valueAt(nRow, nColumn) * m_pB1[nColumn] +
						pIRT->valueAt(nRow, nColumn) * m_pB2[nColumn];
				// J
				m_pSS[centerPt.y() * m_nWidth + rect.left() + nRow].m_cfltIN +=
					pILB->valueAt(nRow, nColumn) * m_pB1[nColumn] +
						pIRB->valueAt(nRow, nColumn) * m_pB2[nColumn];
			}

		// top child
		rectChild.setTopLeft(rect.topLeft());
		rectChild.setBottom(centerPt.y());
		rectChild.setRight(rect.right());
		if (rectChild.width() - 1 != 1 || rectChild.height() - 1 != 1)
			phase2UpdateInnerFlows(rectChild, pLOrT);

		// right child
		rectChild.setTop(rectChild.bottom());
		rectChild.setBottom(rect.bottom());
		if (rectChild.width() - 1 != 1 || rectChild.height() - 1 != 1)
			phase2UpdateInnerFlows(rectChild, pROrB);
	}
}

/****************************************************************************
* End of file												                *
****************************************************************************/

